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o : 

We estimate density of defects frozen into a biological Turing pattern which was 
turned on at a finite rate. A self-locking of gene expression in individual cells, which 
(—^ , makes the Turing transition discontinuous, stabilizes the pattern together with its 

■ defects. A defect-free pattern can be obtained by spatially inhomogeneous activation 

' of the genes. 

OO ' 

Motivation and summary of results 

i— i. 

Long time ago Turing pointed out [1] that simple reaction-diffusion (RD) systems of equations can account for 
formation of biological patterns. The mainstream of research, as reviewed in Ref. [2], is devoted to RD models in 
continuous space. The continuum RD-patterns are smooth and nonpermanent. On the other hand, it is an empirical 
fact that even the nearest neighbor cells can differ sharply in their biological functions and their sets of expressed 
genes. Moreover many biological patterns are permanent. Even the most primitive viruses, like the much studied 
bacteriophage A [3], possess genetic switches that discriminate between different developmental pathways and make 
a once chosen pathway permanent. It is reasonable to assume that cells of higher organisms can also lock their 
£>Y distinctive sets of expressed genes. 

' The Turing patterns on figures in the review [2] are contaminated with defects. If we insist on pattern permanence, 
we must accept that patterns are permanent together with their defects. Sometimes, like for the animal coat patterns, 
permanent defects can provide an animal with its own characteristic life-long but not inheritable "fingerprints". In 
other cases, like formation of vital organ structures, a single defect can be fatal. In this situation it is important to 
!> ■ understand better the origin of defects. 

In this paper we use a simple toy model which in principle should give a homogeneous Turing pattern. Defects 
are particularly manifest on such a simple background. The model has two genes A and B. The genes are strong 
mutual repressors. The strong intracellular mutual inhibition is the factor responsible for pattern permanence. Both 
genes are activated simultaneously in a given cell when a level of their common activator a exceeds its critical value 
a c . Pattern formation in RD models of Ref. [2] was simulated with fixed model parameters. In this paper we turn on 
the activator level a at a finite rate to find a scaling relation between density of defects and the rate. Strong mutual 
intracellular inhibition stabilizes the pattern together with its defects. We obtain permanent domains of A-phase and 
Q ' domains of S-phase divided by sharp cell-size boundaries. 

We also show that an inhomogeneous activation of the genes can result in a perfect defect-free homogeneous pattern. 
At first the activator a exceeds a c in a small seed area where, say, the gene A is chosen. Then a slowly spreads around 
gradually activating more and more cells. The initial choice of A is imposed via intercellular coupling on all newly 
Q-r activated cells. The inhomogeneous activation can be sufficiently characterized by a velocity v with which the critical 
a = a c surface spreads. Thanks to the strong mutual intracellular inhibition there is a nonzero threshold velocity v c , 
such that for v < v c the formation of defects is completely suppressed. In this way the very mutual inhibition which 
is responsible for stability of defects can be harnessed to get rid of them. 

The genetic network that we use in our toy model is functionally equivalent to the genetic toggle switch which was 
syntetized by the authors of the recent paper [4]. In that paper the network is studied experimentally in a single 
"cell". It would be interesting to generalize the experiment to a "multicellular" structure. 
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The toy model 

For the sake of definiteness we take a genetic network with two genes A and B. A and B are mutual repressors. The 
network is symmetric under exchange A «-> B. Expression of both genes is initiated by a common activator a. Let 
A{t, x) and B(t, x) denote time-dependent protein concentrations in the cell at the position x. x belongs to a discrete 
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square lattice with a lattice constant of 1. 
differential equations 



Evolution of the protein concentrations is described by the stochastic 



A(t,x) = RS A (t,x) - A(t,x) 
B(t,x) = RS B (t,x) - B(t,x) 



(1) 
(2) 



The last terms in these equations are responsible for the protein degradation. R is a transcription rate. SA,s(t, x) 6 
{0, 1} are dichotomic stochastic processes. They switch on (0 — > 1) and off (1 — > 0) transcription of a given gene. 
For simplicity the processes are assumed to have the same constant switch-off rate r ofi . The switch-on rates depend 
on concentrations 



r° A n (t,x) = a(t)F 



r° B "(t,x) = a(t)F 



-WB(t,x) + V A ^v) 

n.n.y 

-wA(t,x) + v B (^y) 



(3) 



(4) 



W, V are positive coupling constants, a(t) is a concentration of the activator. F[z] is a smooth step-like sigmoidal 
function; the function F[z] — 10 3 cxp(z — 2.2)/ [1 + exp(z — 2.2)] was used in our numerical simulations. In this model 
the genes A and B are mutual repressors (W > 0). There is a "ferromagnetic" coupling between nearest-neighbor 
cells (V > 0); expression of A in a given cell enhances expression of A in its nearest neighbors. 

The model is motivated by a genetic switch between two mutual repressors like the one studied in the phage A [3] 
and in the E. coli switch [4]. The mutual repressors have a common promoter site on DNA. A necessary condition for 
expression of any of them is a binding of an activator molecule to their promoter site [5] . The concentrations A and 
B influence its affinity to the promoter site. The gene expression is intermittent because of binding and unbinding of 
activator molecules. The nearest-neighbor coupling is possible thanks to signalling through intercellular membrane 
channels. 

In an adiabatic limit, when switching of Sa,b is much faster than protein expression and degradation, the processes 
Sa,b can be replaced by their time averages, 



A(t,x) = 
B(t,x) = 



Ra(t)F 




WB(t,x)+V En.n.yMt,y) 




r oS + a(t)F 


'-WB(t,x)+V En.n.yA&y) 



Ra(t)F 




WA(t,x)+V En.n.yB(t,y) 




r oS + a(t)F 


'-WA(t,x) + V En.n.yB(t,V) 



- A{t,x) . 

- B(t,x) 



(5) 
(6) 



Here we temporarily neglect any noise terms. 



Attractor structure 

In a subspace of uniform configurations A(t),B(t) these equations simplify to the dynamical system 

• _ RaF [-WB + 2dV A] 

~ r off + aF [-W B + 2dV A] ~ ' ( ' 

■ RaF[-WA + 2dV B] _ 



+ aF [-W A + 2dV B] 



where 2d is the number of nearest neighbors in d dimensions. 

The RHS's of these equations define a velocity field on the A — B plane, which is not a gradient field. The 
velocity field has attractor structure which depends on the activator level a. There are two critical activator levels 
a Cl < a C2 . For a < a Cl there is one attractor at [A,B] = [7(a), 7(a)] with an increasing function 7(a). In the range 
a ci < a < a C2 there are three attractors: the old [7(a), 7(a)] plus a new symmetric pair of [a(a), /3(a)] and [/3(a), 01(a)] 
with a(a) > /3(a). For a C2 < a there remain only the two broken symmetry attractors [a(a), (3(a)] and \j3(a),a(a)]. 
The functions a(a),(3(a) and 7(a) are plotted in Fig.l. 
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If we start in the [^4, B] = [0,0] state and slowly increase a- level, the system will stay in the 77-phase until we 
reach a = a C2 . At a = a+ the system will roll into a/3 or /3a-phase. On the other hand, if we start from a C2 < a with 
the system in, say, a/3-phase, then we will have to decrease a down to a = a Cl , where a/3 becomes unstable towards 
the symmetric 77-phase. The discontinuous jumps of the concentrations are illustrated in Fig.l. This hysteresis loop 
is characteristic for first order phase transitions. In the adiabatic limit, where fluctuations are small, there arc no 
short cuts via bubble nucleation. When a Cl (a C2 ) is approached from above (below), the correlation length of small 
fluctuations around this uniform state diverges like in a continuous phase transition. The critical regime is narrow in 
the adiabatic limit so we can rely on the mean field approximation. 



A finite rate Turing transition 



Let us think again about starting from [A,B] = [0,0] and continuously increasing a(t) above a C2 . At a+ the 77 
state becomes unstable and the system has to choose between the a/3 and (3a attractors. If a(t) is increased at a 
finite rate, then there are finite correlated domains which make the choice independently. Despite divergence of the 
correlation length at a~ , the critical slowing down results in a certain finite correlation length £ "frozen" into the 
fluctuations. This scale defines density of defects in the Turing pattern. This effect is well known in cosmology and 
condensed matter physics as Kibble-Zurek scenario [6]. In those contexts the defects disappear rapidly as a result 
of phase ordering kinetics. We will see that in our gene network model the defect pattern is permanent. This effect 
results from a combination of the histeresis loop and the discreteness of the cell lattice. 

To be more quantitative we substitute A(t,x) — 7(a(i)) + SA(t,x) and B(t,x) = j(a(t)) + SB(t,x) into Eqs.(5) 
and linearize them in 5 A, SB. The linearized equations can be diagonalized by (f> — S A — S B and tp = SA + SB. After 
Fourier transformation in space 



/ 



(j>{t,x) = / d d k 4>(t,k) e 



kx 



they become 



<P{t,k)=Rs^(t,k) 



ip(t,k) = Rs^(t,k) + 



r° a Ra(t)F a 
r°ff + a(t)F a ]' 

r oS Ra{t)F' a 
r oS + a(t)F a ] 



W </>(t,k) + V e U (j){t,k) -4>{t,k) , 
-W tp(t, k) + Ve s ip(t, k)] - ip(t, k) 



(9) 

(10) 
(11) 



where = 2^ i=1 cosfci in d dimensions and we skipped the tildas over Fourier transforms. F'[z] — dF[z]/dz and 

we used the shorthands Fa ' = F^>[(—W + 2dV)j(a(t))}. Rs^ are noises which result from fluctuations in RSa,b- 
In the adiabatic limit they can be approximated by white noises (both in space and in time) with small magnitude. 

The next step is to linearize a(t) around its critical value a(t) = a c . 2 + t/r, where r is the transition rate. This 
linearization gives 



r oS Ra(t)F' a t , n2i 

Approximating et — 2d— k 2 in Eqs.(10,ll) and keeping only leading terms in t/r and in k 2 we get 



4>{t,k) = Rs<f,{t,k) + 



C T 



<j>{t, k) 



ip(t,k) = Rs^(t,k) - [2c W + c Vk 2 ] ip{t,k) 



(12) 

(13) 
(14) 



Here we used the identity cq[W + 2dV] = 1, which has to be satisfied because, by definition, 4>(t,0) is a zero mode 
at a C2 . The ip modes are stable for any k. The 0-modes in the neighborhood of k — become unstable for t > (or 
a C2 < a). Eq.(13) is a standard linearized Landau model with the symmetry breaking parameter (ci/co)(t/r) changing 
sign at t = 0. The length scale £ frozen into fluctuations at t > can be estimated following the classic argument 
given by Zurek [6]. For t < the model (13) has an instantaneous relaxation time cor/ci\t\ and an instantaneous 
correlation length c$^JVt / ci\t\. They both diverge at t = 0~. The fluctuations can no longer follow the increasing 
a(t) when their relaxation time becomes equal to the time still remaining to the transition at a = a C21 cor/ci\t\ w \t\. 
At this instant the correlation length is 
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/t/1/2 3/4 \ 

« - (^) • P«) 

This scale determines the typical size of the a/3- and /3a-domains. The scaling relation £ <~ r 1 / 4 was verified by 
numerical simulations illustrated at figures 2 and 3. The domain structures generated in the simulations turned out 
to be permanent. 

The domain structures are permanent because already at a C2 the width of the domain wall interpolating between 
a/3 and /3a is less then the cell size (lattice spacing). The nearest neighbor cells across the wall express different genes. 
The width (the healing length) is determined by the longest length scale of fluctuations around the a/3- or /3a-state. 
These correlation lengths are plotted in Fig. 4. For a > a C2 they are substantially less than 1. In the adiabatic limit, 
where the noises are weak, the domain wall cannot evolve because it would have to overcome a prohibitive potential 
barrier. On a cellular level the barrier originates from the mutual inhibition between A and B in a single cell. Roughly 
speaking, much above a Cl each cell is locked in its gene expression state and insensitive to its nearest neighbors' states. 

Inhomogeneous activation 

The intracellular mutual inhibition stabilizes the Turing pattern but it also stabilizes the defects frozen into the 
pattern. With the £ ~ t 1 / 4 scaling the number of defects is rather weakly dependent on r. There may be not enough 
time during morphogenesis to get rid of the defects by simply increasing r. However, it is possible to generate a 
defect-free pattern by spatially inhomogeneous switching of the activator level a. For example, its concentration can 
exceed a C2 at one point at first, where the cells happen to pick (or are forced to pick), say, a/3-phase, and then the 
activator can gradually spread around so that the initial seed of a/3-cells gradually imposes their choice on the whole 
system. For continuous transitions this effect was described in Ref.( [7]). 

The effect of defect suppression in inhomogeneous activation can be most easily studied in a one dimensional version 
of the model (1). Suppose that a smooth activator front is moving across the one dimensional chain of cells with a 
velocity v, a(t,x) « a C2 + (vt — x)/vt close to x = vt where a = a C2 . For definiteness we impose two asymptotic 
conditions: for vt <C x ( where a < a C2 ) the cells are in the 77-state, and for x <C vt (where a > a C2 ) they are in 
a/3-phase. We can expect that as the a-front moves to the right it is followed by the a/3 front gradually entering the 
area formerly occupied by the 77-phase. If the concentration front is fast enough to move in step with the activator 
front, then the a/3-phase will gradually fill the whole system. If, on the other hand, the concentration front is slower 
than the activator front then the front of the a/3-phase will lag behind the a = a C2 front. The gap between the two 
fronts will grow with time. The gap will be filled with the unstable 77-phase (a > a C2 behind the a-front). When 
the gap becomes wide enough, then 77-state will be able to decay towards the /3a-state. A domain of /3a-phasc 
will eventually be nucleated behind the a-front. Now the /3a-domain will grow behind the a-front until its front lags 
sufficiently behind so that a new domain of a/3-phase will be nucleated. In this way the activator front will leave 
behind a landscape of alternating a/3- and /3a-domains qualitatively the same as for homogeneous activation. 

The success of the inhomogeneous activation depends on the relation between the velocity v of the a-front and that 
of the concentration front. As illustrated in Fig. 4 fluctuations around the a/3-state have two families of modes each 
with a different correlation length. For any a each fc-mode within each family has a different diffusion velocity: a ratio 
of its wavelength to its relaxation time. The lowest of these diffusion velocities, v c (a), is the maximal velocity at which 
the a/3-phase can spread into the area occupied by the 77-phase. v c (a C2 ) =v C2 > because at a = a C2 the a/3-state is 
stable (the hysteresis loop again!). v c (a) increases with an increasing a. If v < v C2 the a/3-front moves in step with the 
a-front; its tail spreads into the vt < x area imposing an a/3-bias on the fluctuations around 77-state. The a/3-phasc 
spreads without nucleation of any /3a-domains. For v < v C2 a defect-free uniform Turing pattern forms behind the 
activator front. Results from numerical simulations of the inhomogeneous activation are presented in Fig. 5. 

More complicated patterns 

Finally, it is time to comment on more complicated models which are expected to give more complicated patterns 
than the (in principle) uniform pattern discussed so far. Let us pick a zebra pattern for example. For the uniform 

pattern the first mode to become unstable in Eq.(13) is the k = mode. The final pattern has an admixture of 
k's in a range w around k = 0. In distinction, for the zebra pattern the first unstable modes are those on the 
circle \k\ = 2n/L, where L is the spacing between zebra stripes. The final pattern has an admixture of k's in a ring 
of thickness « around the circle \k\ = 2tt/L, compare results for Swift- Hohenberg equation in Ref. [8]. This 
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admixture results in defects frozen into zebra pattern. The inhomogcneous activation can be applied in the zebra case 
too. In addition it can be used to arrange the stripes. An activator spreading from an initial point would result (at 
least close to the initial point) in concentric black and white rings. A front of activator moving through the system 
would comb the stripes perpendicular to the front. 
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FIG. 1. The thick lines are: a(a) (top), 7(a) (middle), /9(a) (bottom). The vertical lines with arrows illustrate the discon- 
tinuous jumps by the concentrations A and B during the a/3 — > 77 transition at o cl « 0.34, and the 77 — ► a/3 transition at 
a C2 w 0.37. Model parameters used in this graph are: R = 4, W = 3, V = 1, d = 2, r° a = 10 3 . 
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FIG. 2. Permanent pattern obtained after switching-on the activator on a 32 x 32 periodic lattice. It is a contour plot of 
A — B; white is A-rich (af3) and black is B-rich (Pa). The activator was turned on as a(t) = t/r with r = 32 and t € (0, 32). 
Model parameters were the same as in Fig.l. A discrete time step was At = 1CF 4 . 
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FIG. 3. log(£) as a function of log(r). £ was obtained as an average domain size along a cross section through patterns like 
that in Fig. 2. For any given r the average was taken over outcomes of many simulations and over all the possible vertical and 
horizontal cross sections. The vertical point size is a triple standard deviation. The simulations were done on a 1024 x 1024 
lattice. The slope was fitted as 0.24 ± 0.02, which is consistent with the predicted 0.25. 
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FIG. 4. The correlation lengths of the fluctuations around the state a/3 as functions of a. The vertical gridlines mark 
a ci ~ 0.332 and a C2 » 0.375. The larger correlation length diverges at a ci . These correlation lengths should be compared 
with the lattice spacing which is 1. The correlation lengths were obtained by expanding A(t,x) = a(a) + 8A(t,x) and 
B(t,x) = /3(a) + SB(t,x), Fourier-transforming the fluctuations in space and subsequent diagonalization for small k. 
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FIG. 5. Density n of domain walls between a/3 and /3a-states behind an activator front with velocity v. The activator was 
a(t,x) = (vt — x)/vt for x < vt and a(t,x) = for vt < x. 1/vt — 0.1 was kept fixed so that the slope of a versus x was 
independent of v. The model parameters were the same as in Figs. 1,2 but with d = 1 instead of 2 and V = 2 instead of 1 
{Vd = 2 as before). For these model parameters v c w 0.9 in consistency with the numerical results. 
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